Estimating the benefit of green infrastructure to urban ecosystems: A synthesis and case-study

Alessandro Filazzola

Scott MacIvor - UTSC

Namrata Shrestha - TRCA

Glenn Milner - OCC

library(tidyverse)
library(PRISMAstatement)

Global urbanization continues unabated, with more than 50% of the worlds’ population living in cities. Cities are conventionally viewed as a threat to local biodiversity because natural habitat is replaced with development. However, more recently, there is greater acknowledgement from the public and private sectors that supporting local environments sustains critical ecosystem services, which in turn improves human health and biodiversity conservation. Consequently, urban planning and design has shifted towards green infrastructure (GI), such as green roofs and retention ponds, to increase connections between city and nature in an era of climate change. The contribution of GI to some ecosystem services has been proven (e.g. stormwater management, building cooling), but the contribution to biodiversity conservation remains unspecified. Using a systematic literature review, this project will (i) determine effect estimates that relate different GI types and characteristics to the impacts on natural systems, and (ii) compile relevant data to develop different implementation scenarios GI for Toronto and region. This study will inform natural system planning and improve quantification of GI on urban ecosystems. Findings from this research will have global ramifications that allow city planners to optimize GI implementation for sustainable development and decrease the impacts of cities on natural systems.

Objectives

1. A meta-analysis of the literature summarizing the effect of green infrastructure on natural systems.
2. Using Toronto, Ontario as a case study, develop a tool that can communicates the effectiveness of different GI implementation for biodiversity conservation.

Expected Deliverables

1. A peer-reviewed journal article that is a meta-analysis from objective 1.
2. A tool or data analysis that projects different scenarios of green infrastructure implementation for the City of Toronto.

Timeline

date	task
June 18	Begin meeting with staff and MacIvor lab to set out workplan
June 25	Begin literature review and data extraction
July 2	Aggregate available data for GI analysis in Toronto
July 3	Complete meetings with TRCA staff on relevant considerations for the project
July 9	Determine important parameters for modelling GI in Toronto
August 20	Complete collection and review of relevant articles
September 3	Conduct meta-analysis on available data
September 10	Propose candidate models for quantifying GI effects for natural systems
Sept 24	Model validation and begin writting manuscript
Oct 15	Complete a draft of manuscript and finalize model

Literature Review - 1. Search

A systematic literature search will be conducted using Web of Scienc for all articles between 1980 and 2018. This time frame was chosen because it captures the majority of the literature on green infrastructure and with articles that have some measured estimates. The review will include all studies globally. The intended purpose of this search is to capture all articles that have documented both green infrastructure implementation and a measure of natural systems. The search terms that will be used are:

*green infrastructure OR low*impact development OR Sustainable Drainage System* OR Water Sensitive Urban Design OR green*roof AND *diversity OR species OR ecosystem OR ecology OR habitat* OR co-benefit

These terms have returned 871 results (as of June 2018)

Revise list

search1.1 <- read.csv("data/WOS-lit.csv")
search1.2 <- read.csv("data/WoSPart3-July_4_2018.csv")
net.difference <- anti_join(search1.2, search1.1, by = "DOI")
net.difference <- net.difference %>% select(Title, DOI) #to simplify for a look
nrow(net.difference) #count of number of differences from consecutive search

## [1] 182

## 182 papers to be added by including revised terms
## Select those articles and join with other dataset
net.difference <- anti_join(search1.2, search1.1, by = "DOI")

updated.search <- rbind(search1.1, net.difference)

#write.csv(updated.search, "data/WOS-lit.updated.csv")

Adding revised terms from July 3rd meeting added 182 papers Total articles returned = 1,053 (as of July 2018)

## Adding terms for naturalized pond and pollinator garden
search1.2 <- read.csv("data/WOS-lit.updated.csv")
search1.3 <- read.csv("data/WoSPart4-July_11_2018.csv")
net.difference <- anti_join(search1.3, search1.2, by = "DOI")
net.difference <- net.difference %>% select(Title, DOI) #to simplify for a look
nrow(net.difference) #count of number of differences from consecutive search

## [1] 0

## 213 papers to be added by including revised terms
## Select those articles and join with other dataset
net.difference <- anti_join(search1.3, search1.2, by = "DOI")

updated.search <- rbind(search1.2, net.difference)

#write.csv(updated.search, "data/WOS-lit.updated.csv")

Adding revised terms from July 3rd meeting added 213 papers Total articles returned = 1,224 (as of July 2018)

Literature Review - 2. Sort

This steps includes a. checking for duplicating, b. reviewing each instance for relevancy, c. consistently identifying and documenting exclusion criteria. Outcomes include a list of publications to be used for synthesis, a library of pdfs, and a PRISMA report to ensure the worflow is transparent and reproducible. Papers were excluded with the following characteristics:

• Not emperical study (e.g. review, book chapter)
• Irrelevant categories (e.g. political science, law, sports tourism, art)

evidence <- read.csv("data/evidence.csv")

excludes <- evidence %>% group_by(reason) %>% count(exclude) %>% filter(reason!="")
ggplot(excludes, aes(x=reason, y=n)) + geom_bar(stat="identity") + coord_flip()

## frequency of study
year.rate <- evidence %>% group_by(Year) %>% summarize(n=length(Year))

Papers processed - Progress

## Completed so far
prog <- sum(evidence$exclude!="")
prog

## [1] 1266

## Remaining
total <- nrow(evidence)
total

## [1] 1266

setTxtProgressBar(txtProgressBar(0,total,  style = 3), prog)

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

Initial pass for relevant papers complete.

Description of studies

GI.type <- evidence %>% group_by(GI.type) %>% count(exclude) %>% filter(GI.type!="")
ggplot(GI.type, aes(x=GI.type, y=n)) + geom_bar(stat="identity") + coord_flip()

Representations of relevant GI types found in papers

Prisma report

## total number of papers found
nrow(evidence)

## [1] 1266

## number of papers found outside of WoS
other <- read.csv("data/other.sources.csv")
nrow(other)

## [1] 28

## number of articles excluded
excludes <- evidence %>% filter(exclude=="y")
nrow(excludes)

## [1] 1123

## relevant papers
review <- evidence %>% filter(exclude!="y")
nrow(review)

## [1] 143

## papers for meta
meta <- evidence %>% filter(meta.=="yes")
nrow(meta)

## [1] 117

prisma(found = 1255,
       found_other = 28,
       no_dupes = 1283,
       screened = 1283,
       screen_exclusions = 1140,
       full_text = 143,
       full_text_exclusions = 0,
       qualitative = 143, 
       quantitative = 66,
       width = 800, height = 800)

## Loading required namespace: DiagrammeR

Literature Review - 3. Synthesis

The research questions we are exploring:

1. What are the patterns of GI studies globally
2. How does green infrastructure compare to conventional “grey” equivalents (e.g. green roof to conventional roof)?
3. How does green infrastructure compare to its natural equivalents (e.g. retention ponds )?
4. What features of green infrastructure can improve the quality of natural systems?

Patterns of GI Studies Globally

require(ggmap)
###  Start with base map of world
mp <- NULL
mapWorld <- borders("world", colour="gray50", fill="gray50") # create a layer of borders
mp <- ggplot() +   mapWorld

## colorblind-friendly palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#000000")

meta <- read.csv("data//evidence.csv")
meta <- subset(meta, GI.type!="")

## plot points on top
mp <- mp+ geom_point(data=meta , aes(x=lon, y=lat, color=GI.type), size=2) + scale_colour_manual(values=cbPalette)+
    theme(legend.position="bottom", text = element_text(size=20))
mp

## Number of studies extracted from online data
occurdat<-list.files("data//MS.data",pattern=".csv$",full=T)
length(occurdat)

## [1] 76

## 70 Studies found with usable data for synthesis

Green infrastructure comparison to conventional

## load master datasets
meta <- read.csv("data//Master.GI.Datasets.csv")
## Omit repo 3 because duplicated with study 1305 and remove repo-9 because not equivalent GI comparisons. Removed Repo 3 because compare roof with ground
meta <- meta %>% filter(Study != "repo-3" & Study!="repo-9" & Study!="repo-1") 

## Drop relative abundance because difference = 0 
meta <- meta %>% filter(Estimate!="Relative.Abundance")

## Load packages and functions
library(reshape2)
library(metafor)
source("meta.eval.r") ## Multiple aggregate


## Create Unique identifier column
meta2 <- meta
meta2[,"UniqueSite"] <- paste(meta2$Study, meta2$Taxa.simplified, meta2$GI.compare, meta$Estimate, sep="-")
meta3 <-  meta2 %>% filter(Infrastructure != "natural") %>%  filter()

## Use function to extract summary statistics for comparisons
## meta.eval  arguments are (meta.data, compare, ids , stats)
Infra.compare <- meta.eval(meta3, Infrastructure, UniqueSite, Stat )

## Combine the lists into same dataframe
## Rename Columns in second dataframe
Infra.stat <- Infra.compare[[2]] ## extracted statistics 
names(Infra.stat) <- c("UniqueSite","green_mean","green_sd","grey_mean","grey_sd","green_n","grey_n") ## rename columns to match
Infra.raw <- Infra.compare[[1]] ## calculated statistics from raw values

## Join two dataframes
meta.stat <- rbind(Infra.raw, Infra.stat[, names(Infra.raw)])


meta.ready <- escalc(n1i = grey_n, n2i = green_n, m1i = grey_mean, m2i = green_mean, sd1i = grey_sd, sd2i = green_sd, data = meta.stat, measure = "SMD", append = TRUE)

taxa.test <- c(rep("invertebrate",5),"amphibian","birds.bats","birds.bats","invertebrate","invertebrate","birds.bats","birds.bats","invertebrate")
GI.test <- c("wall","wall","roof","roof","pond","pond","wall",rep("roof",5),"roadside")


#random-effects meta-analysis for green infrastructure vs grey
m1 <- rma(yi=yi, vi=vi, mods=~factor(taxa.test), data = meta.ready)
summary(m1) 

## 
## Mixed-Effects Model (k = 13; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -26.8388   53.6776   61.6776   62.8879   69.6776  
## 
## tau^2 (estimated amount of residual heterogeneity):     10.9891 (SE = 5.0836)
## tau (square root of estimated tau^2 value):             3.3150
## I^2 (residual heterogeneity / unaccounted variability): 98.11%
## H^2 (unaccounted variability / sampling variability):   52.99
## R^2 (amount of heterogeneity accounted for):            11.12%
## 
## Test for Residual Heterogeneity: 
## QE(df = 10) = 128.5068, p-val < .0001
## 
## Test of Moderators (coefficient(s) 2:3): 
## QM(df = 2) = 3.8058, p-val = 0.1491
## 
## Model Results:
## 
##                                estimate      se     zval    pval     ci.lb
## intrcpt                         -1.7060  3.3264  -0.5129  0.6080   -8.2256
## factor(taxa.test)birds.bats     -3.3826  3.7411  -0.9042  0.3659  -10.7150
## factor(taxa.test)invertebrate    0.6600  3.5310   0.1869  0.8517   -6.2607
##                                 ci.ub   
## intrcpt                        4.8136   
## factor(taxa.test)birds.bats    3.9498   
## factor(taxa.test)invertebrate  7.5806   
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Produce a forest plot to determine the effect sizes for each study
forest(m1, slab = meta.stat$UniqueSite)

## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m1)

Green infrastructure comparison to natural equivalent

## Create Unique identifier column
meta2 <- meta
meta2[,"UniqueSite"] <- paste(meta2$Study,  meta2$Taxa.simplified, meta2$Nat.compare, meta2$Estimate, sep="-")

## Remove comparisons except urban and rural
meta2 <- meta2 %>% filter(Habitat == "urban" | Habitat == "natural")

## Determine the number of comparisons available 
compare.eval(meta2, Habitat, UniqueSite)

## [1] 20

## Use function to extract summary statistics for comparisons
## meta.eval  arguments are (meta.data, compare, ids , stats)
nat.compare <- meta.eval(meta2, Habitat, UniqueSite, Stat )


## Combine the lists into same dataframe
## Rename Columns in second dataframe
nat.stat <- nat.compare[[2]] ## extracted statistics 
names(nat.stat) <- c("UniqueSite","natural_mean","natural_sd","urban_mean","urban_sd","natural_n","urban_n") ## rename columns to match
nat.raw <- nat.compare[[1]] ## calculated statistics from raw values

## Join two dataframes
meta.stat <- rbind(nat.raw, nat.stat[, names(nat.raw)])

meta.ready <- escalc(n1i = urban_n, n2i = natural_n, m1i = urban_mean, m2i = natural_mean, sd1i = urban_sd, sd2i = natural_sd, data = meta.stat, measure = "SMD", append = TRUE)

#random-effects meta-analysis for urban GI vs natural
m1 <- rma(yi, vi, data = meta.ready)
summary(m1) 

## 
## Random-Effects Model (k = 20; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -14.8475   29.6950   33.6950   35.5839   34.4450  
## 
## tau^2 (estimated amount of total heterogeneity): 0.0758 (SE = 0.0539)
## tau (square root of estimated tau^2 value):      0.2753
## I^2 (total heterogeneity / total variability):   49.98%
## H^2 (total variability / sampling variability):  2.00
## 
## Test for Heterogeneity: 
## Q(df = 19) = 39.3013, p-val = 0.0040
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb   ci.ub   
##  -0.0889  0.0957  -0.9295  0.3526  -0.2764  0.0986   
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Natural vs Urban GI
## Produce a forest plot to determine the effect sizes for each study
forest(m1, slab = meta.stat$UniqueSite)

## Check for publication bias
## The symetrical distriubtion suggests there is no publication bias
funnel(m1)

Green infrastructure features in relation to measures of green infrastructure

## The area of green infrastructure
names(meta)[20:21] <- c("GI.area","height")
meta.area <- subset(meta, GI.area>0)

## Determine unique identifier
meta.area[,"UniqueSite"] <- paste(meta.area$Study, meta.area$Taxa.simplified, meta.area$Nat.compare, meta.area$Estimate, sep="-")

## Summarize average richness by area sizes
area.stat <- meta.area %>% filter(Stat=="count" | Stat=="mean") %>% filter(Estimate=="richness") %>%   group_by(GI.area, GI.type) %>% summarize(val=mean(Value)) %>%  data.frame(.)

area.rich <- area.stat %>% filter(GI.type=="green roof" | GI.type=="green wall" | GI.type=="retention pond" | GI.type=="yards/home gardens" | GI.type=="public/community gardens")

library(ggplot2)

## Species richness per area
ggplot(area.rich,  aes(x=GI.area, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Species Richness") + xlab("Average Area of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## Summarize average abundance by area sizes
area.stat <- meta.area %>% filter(Stat=="count" | Stat=="mean") %>% filter(Estimate=="abundance") %>%   group_by(GI.area, GI.type) %>% summarize(val=mean(Value)) %>%  data.frame(.)

area.abd <- area.stat %>% filter(GI.type=="green roof" | GI.type=="green wall" | GI.type=="retention pond" | GI.type=="yards/home gardens" | GI.type=="public/community gardens") %>% 
  filter(GI.area<50000) ## keep numbers approximately similar - removed outlier of 200,000

## Species richness per area
ggplot(area.abd,  aes(x=GI.area, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Abundance") + xlab("Average Area of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## Compare pH for abundance
ph.data <- meta %>%  filter(pH>0) %>%  filter(GI.type=="retention pond" | GI.type=="natural water")
ph.abd <- ph.data %>% filter(Estimate=="abundance" & Stat=="count")%>%   group_by(pH, GI.type) %>% summarize(val=mean(Value)) %>%  data.frame(.)

## plot richness against pH
ggplot(ph.abd,  aes(x=pH, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Richness") + xlab("Average pH of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## Compare pH for occurrence
ph.occ <- ph.data %>% filter(Estimate=="occurrence")%>%   group_by(pH, GI.type) %>% summarize(val=mean(Value)) %>%  data.frame(.)

m1 <- glm(Value~ pH, family=binomial, data=ph.data %>% filter(Estimate=="occurrence"))
summary(m1)

## 
## Call:
## glm(formula = Value ~ pH, family = binomial, data = ph.data %>% 
##     filter(Estimate == "occurrence"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8915  -0.8488  -0.8377   1.5229   1.6131  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.22413    1.45261  -0.154    0.877
## pH          -0.08643    0.20324  -0.425    0.671
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 308.55  on 251  degrees of freedom
## Residual deviance: 308.37  on 250  degrees of freedom
## AIC: 312.37
## 
## Number of Fisher Scoring iterations: 4

## plot richness against pH
ggplot(ph.occ,  aes(x=pH, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Occurrence") + xlab("Average pH of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## plot richness against height
h.rich <- meta %>%  filter(height>0)  %>% filter(Estimate=="richness") %>% filter(Stat=="count" | Stat=="mean") %>% group_by(height, GI.type) %>%  summarize(val=mean(Value)) %>%  data.frame(.)

ggplot(h.rich,  aes(x=height, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Richness") + xlab("Average Height of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## plot abundance against height
h.abd <- meta %>%  filter(height>0)  %>% filter(Estimate=="abundance") %>% filter(Stat=="count" | Stat=="mean") %>% group_by(height, GI.type) %>%  summarize(val=mean(Value)) %>%  data.frame(.)

ggplot(h.abd,  aes(x=height, y=val, color=GI.type) ) + geom_point(size=3) + theme_bw() + scale_color_brewer(palette="Set2") + ylab("Average Abundance") + xlab("Average Height of Green Infrastructure")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

## Pond salinity
salt <- subset(meta, Salinity>0)
length(unique(salt$Salinity)) ## not enough samples for a meaningful comparison

## [1] 4

## Pond depth
deep <- subset(meta, depth..m.>0)
length(unique(deep$depth..m.)) ## not enough samples for a meaningful comparison

## [1] 4

Case Study - 4. Green infrastructure Mapping

gf.data <- read.csv("data//GI.data//GreenRoofGeocoded.csv")
garden.data <- read.csv("data//GI.data//ComGardensdata.csv")
pond.data <- read.csv("data//GI.data//RetentionPondsGPS.csv")

## Extract coordinates only for GI
gf.data <- gf.data[,c("lon","lat")]
garden.data <- garden.data[,c("lon","lat")]
pond.data <- pond.data[,c("lon","lat")]

## combine into single dataset
GI.data <- rbind(gf.data,garden.data, pond.data)
## add column for data type
GI.data[,"GI.type"] <- c(rep("Green Roof",nrow(gf.data)),rep("Community Garden",nrow(garden.data)),rep("Retention Pond", nrow(pond.data)))

## interactive map
library(leaflet)

## Warning: package 'leaflet' was built under R version 3.4.4

m <- leaflet(data=GI.data) %>%
  addTiles() %>%  
   addCircleMarkers(~lon, ~lat, color =  ifelse(GI.data$GI.type == "Retention Pond", 'blue', 'red'), radius=5)
m